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When classical systems fail to explore their entire configurational space, intriguing macroscopic 
phenomena like aging and glass formation may emerge. Also closed quanto-mechanical systems 
may stop wandering freely around the whole Hilbert space, even if they are initially prepared into 
^ , a macroscopically large combination of eigenstates. Here, we report numerical evidences that the 

D ■ dynamics of strongly interacting lattice bosons driven sufRciently far from equilibrium can be trapped 

into extremely long-lived inhomogeneous metastable states. The slowing down of incoherent density 
\ excitations above a threshold energy, much reminiscent of a dynamical arrest on the verge of a glass 

■ transition, is identified as the key feature of this phenomenon. We argue that the resulting long-lived 
I I' inhomogeneities are responsible for the lack of thermalization observed in large systems. Such a rich 

' phenomenology could be experimentally uncovered upon probing the out-of-equilibrium dynamics 

of conveniently prepared quantum states of trapped cold atoms which we hereby suggest. 

a: 

I _ The ergodicity axiom in classical statistical mechanics states that, during its time evolution, a closed macroscopic 

■ system uniformly explores the entire phase space compatible with conservation laws, so that the time average of 
\ any observable comes to coincide with the micro-canonical ensemble average and, when the observable is local, also 

I t ■ with the canonical Gibbs ensemble average. Nonetheless, ergodicity can be violated in classical systems, a noticeable 
\ example being glasses. [l[ Quantum effects might also spoil ergodicity by preventing the wave function from diffusing 
a ' within all available configurations. This phenomenon is actually known to occur in the presence of disorder and 
manifests itself either by single-particle Q or many-particle [^-Q wave function localization. However, alike classical 
^ models for glassy behavior, ergodicity breakdown in the quantum dynamics may not necessarily require disorder and 
Q \ it could instead be entirely due to frustrating dynamical constraints. 0, 0] This issue is currently attracting great 
O ■ interest, since well controlled realizations of closed quantum systems have become feasible upon trapping cold 

I atomic species. [lo| Indeed, similarly to what can be done in numerical simulations, one can prepare atoms in a 
^vq i given initial state and probe their time evolution under a Hamiltonian whose parameters are fully under control, thus 



^ ] offering the unique opportunity to monitor the ergodicity principle at work in the quantum realm. [11 

. In this work, we report numerical evidences that an isolated system of strongly interacting bosons, modeling atoms 
' ' ■ in optical lattices, can be trapped during its evolution into long-lived inhomogeneous metastable states, provided that 
. its internal energy exceeds a certain threshold. We argue that the slowing down of high-energy incoherent excitations 
in the strongly correlated system is the key feature responsible for this dynamical arrest, much resembling a kind of 
glass transition. By formulating the problem in a different language, we explicitly show that a system initially prepared 
in a inhomogeneous state is unable to diffuse within the entire configurational space; such a dynamical localization in 
the many-body Hilbert space looks intriguing and may represent a kind of many-body Anderson localization that 
occurs without disorder. The above phenomenon is put in further relation with and deemed responsible for the lack of 
ergodicity observed in large finite size systems. This belief is confirmed by means of a novel time-dependent variational 
Monte Carlo method that we introduce hereby and an experimental set up to uncover this very rich phenomenology 
is also suggested. 

One of the simplest models that can be realized in experiments is the Bose-Hubbard Hamiltonian: l^, 12 1 

characterized by the amplitude J for a bosonic species of an atom to hop between nearest-neighboring wells of an 
optical lattice and by a local repulsion U among atoms localized in the same potential well. The operators b\ and 
create and destroy, respectively, a boson on site i, and Ui ~ b\b^ is the density operator, [s^l Experiments are often 
performed with anisotropic lattices that realize a collection of almost uncoupled chains, a fortunate case for numerical 
simulations that we shall mainly consider hereafter, apart from a brief excursion in two dimensions towards the end 
of the paper. 
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Figure 1: Inhomogeneous initial states constituted by clusters of doubly occupied sites (large blue circles) and empty sites 
(small dots). The size of each cluster is denoted by Sd- 



In one dimension, obstacles to ergodicity can arise in intcgrable models. [13[ However, the Hamiltonian ([T} is not 
integrable and, indeed, there are experimental evidences that its dynamical evolution may succeed in fast relaxing to a 
thermal state. Specifically, in a recent experiment a system of ®''Rb atoms, well described by the Hamiltonian ff]) , has 
been prepared in a state in which the sites of the optical lattice were alternatively empty and singly occupied. |14| This 
state was let evolve for different experimental conditions corresponding to different ratios U/ J. Even at the largest 
value U/ J ~ 10, the initial density profile (...1,0,1,0,...) was found to rapidly relax to the homogeneous thermal 
one (. ..^,^,^,^,...), with half a boson per site, much faster than the integrable counterparts of non- interacting or 
infinitely- interacting (i.e., hard-core) bosons and consistently with the increased number of relaxation channels that 
opens once integrability is lost. 

Numerical simulations of the above experiment successfully reproduce the observed thermal behavior, [l^l However, 
there are also numerical evidences pointing out a breakdown of ergodicity in the same model ([T]) but within a different 
region of the parameter space, specifically when the number of bosons is one per site. This case, as well as any 
other one with integer density, is special because, at equilibrium and at zero temperature, the model ([T]) undergoes a 
quantum phase transition into a Mott insulator above a critical U/ J. Even though the phase transition is washed out 
by thermal fluctuations, nonetheless its influence on the spectrum seems to prevent thermalization above a certain 
U/J in numerical simulations of finite size systems[2^ [2^; a result that may or may not [2^ prelude to a true 
breakdown of ergodicity in the actual thermodynamic limit. 

Here, we propose an experiment a lot alike the one previously described, (3] which may distinguish very sharply a 
change from chaotic to non-chaotic dynamics in the Bose- Hubbard model ([T]). Our proposal starts by observing that, 
in the gapless phase next to the Mott insulator at one boson per site, low-energy itinerant Bogoliubov quasiparticles 
should coexist with high-energy incoherent excitations, which, for U/J 3> 1, can be identified as sites occupied by 
more than a single boson. It is well possible that the overall relaxation depends critically upon the population of 
those high-energy excitations, [l^l which can be assessed by tailoring initial states with lots of doubly occupied sites 
rather than none as in the experiment of Ref . [l^ . This is indeed what we propose and simulate by means of different 
and complementary numerical tools, such as exact diagonalization, time-evolving block decimation [isj and a novel 
time-dependent variational Monte Carlo algorithm, as discussed in more detail in the Supplementary Material. 



Results 



Inhomogeneous initial states and dynamical localization 

We start from analyzing the model at density n — 1, when at equilibrium a Mott transition occurs at a critical 
(U/J)c — 3.5. [3| We imagine to prepare an initial state where all the sites are either empty or doubly occupied. In 
particular we consider the states depicted in Fig. [TJ namely with clusters of doubly occupied sites of variable size Sd- 
These states are let evolve with the spatially homogeneous Hamiltonian ([1]) for different U/J, below and above the 
critical value. While for small U/J the density profile rapidly reaches the equilibrium configuration (...1,1,1,1,...), 
for large U/J, it stays close to its initial value for a remarkably long time. Eventually, since the system is finite, the 
density profile approaches the homogeneous plateau, with small residual oscillations that get damped as the system 
size increases. 

We can define a relaxation time tr (whose inverse is shown in Fig. [2l for Sd — 1), as the first time for which the 
local density approaches its homogeneous value. We highlight that tr, at a specific ([//J)^^", has a sudden step up, 
which becomes sharper and sharper as the system size increases. The above results show that above (C//J)^^" the 
system has the tendency to stay dynamically trapped into long-lived inhomogeneous configurations. 
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Figure 2: Left panel - Inverse relaxation times of the local density for the initial state (. . . 2, 0, 2, 0, . . . ). Exact diagonal- 
ization results are reported, with darker points marking larger systems, respectively A*" = 8, 10 and 12 with periodic-boundary 
conditions. In the Inset, the time dependence of the on-site densities is shown. Right panel - Effective-Hamiltonian evolution 
of the average density within a cluster of different size Sd, where the dimensionless time is defined as t''^^ = ^^t. We show 
results for N = 72 and 96 sites with open-boundary conditions, obtained by the time-evolving block decimation technique. [l3| 
We note that even small clusters of doubly occupied sites can effectively freeze the dynamical evolution in the large-C/ regime. 



We emphasize that such a behavior is all the more remarkable not because doubly occupied sites seem unable to 



decay, which is known to require for U / J ^ 1 very long times and large system sizes, [18| but rather because they 
are not capable to move, hence restore translational symmetry. We can better understand this surprising result by 
an effective Hamiltonian that can be derived following the same reasoning of Refs. [l^ H^- For a sufficiently large 
interaction, it is justified to project the evolution onto states with the same potential energy per site C//2, at least for 
time scales shorter than / . One reafizes that states with the same potential energy but with triply and singly 
occupied sites start to contribute only at order J^'/C/^, so that, with accuracy J^/i7, the number of doubly occupied 
sites is conserved. If we associate a fictitious spin up or down to a doublon (doubly occupied site) or a holon (empty 
site), respectively, we find that the effective Hamiltonian that controls the evolution reads: 

(2) 



which describes a hard-axis ferromagnetic Heisenberg model. [19[ We note that the evolution of an XXZ spin chain 
starting from an antiferromagnetic ordered initial state has been studied numerically in Ref. plj with DMRG. Here 
we consider the effective dynamics ([2]) in a regime when the anisotropy would favor a ferromagnetic ground state and 
when the initial state contains clusters of up and down spins of increasing lenght. The smaller cluster size of 5^ = 1 
corresponds to a Neel initial state considered in Ref [2lj. In the right panel of Fig. [5] we show the time-dependence 
of the clusters density in the large U regime, evolved according to the effective Hamiltonian ©. Remarkably, we see 
that even for small clusters the system fails to restore the spatial homogeneity up to very long time scales, which turn 
to be far beyond those currently accessible in typical experimental setups. 

The slowing down of the dynamics can be traced back to the effective attraction among doublons, [2^ which makes 
their aggregates hard to break up. In other words, what seems to matter more is the dissociation of clusters of 
doublons, rather than the decay of a single one. Indeed we have checked (not shown) that in the large U regime 
the system has the tendency to get dynamically stuck into clusters of doublons of finite size, whereas a much faster 
annihilation and recombination rate is observed in the small U limit. 

To get further insights into the dynamical behavior of the system, we recast the problem in a different language. 
Starting from the initial state, denoted as [0), we can generate an orthogonal basis set i = 0, 1, . . . , by repeatedly 
applying the Hamiltonian (see Supplementary Material). In this Lanczos basis of many-body wave functions, the 
Hamiltonian has the form of a tight-binding model on a semi- infinite chain. Each site i = 0, 1, . . . corresponds to 
a many-body state, it has an on-site energy = and is coupled only to its nearest neighbors by hopping 

elements ti^i+i and ti^i-i. (22j It is easy to realize that the unitary evolution of the original many body problem is 
thus fully equivalent to the dynamics of a single particle, initially sitting at site 0, that is then let propagate along 
such a tight-binding chain of many-body states. We note that, both ti and U^i+i largely fiuctuate from site to 
site, therefore resembling an effective Anderson model, even though those parameters are in reality deterministic, see 
Fig. [3] In the same figure, we also show the mean distance traveled by the particle after time t starting from the 



4 




Figure 3: Left panels - On site energies and nearest-neighbor hoppings of the effective chain that represents the Hamiltonian in 
the Lanczos basis starting from the state (. . . , 2, 0, 2, 0, . . . ). Red points refer to U — 2 J, wlien tire particle does diffuse starting 
from site 0, while blue points to U = lOJ, when it does not. The shaded regions correspond to energies less or equal that of 
the initial state. Right panel - Time-dependent expectation value of the wave-packet position of the effective particle traveling 
in the Hilbert space generated by a chain of Timax = 1000 Lanczos states. The red points correspond to U = 2J and the blue 
points to U = 10 J (the original lattice size of the Bose-Hubbard model is = 12). The shaded region marks the center of the 
Lanczos chain, which is not reached in the localized regime. 



first site of the chain, which corresponds to an initial state |0) = (. . . 2, 0, 2, 0, . . . ), and for different U/ J. We observe 
that, for small values of J7/ J, the particle diffuses and its wave-packet finally spreads over the whole chain, in a rather 
uniform way. On the contrary, above a certain critical value of the interaction ([//J)J?^", the particle stays localized 
near the origin for arbitrarily long times. A closer look to the structure of the on site energies reveal the existence of 
a potential well at the edge of the chain. This is crucial in order to understand the observed localization transition. 
Since the potential well cannot induce a true bound state below the bottom of the spectrum, at most a resonance 
will form in the spectrum, from which the particle could in principle escape in a finite time. This indeed happens 
at small U but apparently not at large U, where the increased depth of the well and the effective randomness of the 
on-site energies conspire together to keep the particle localized close to the edge, preventing the the states in the well 
to hybridize with other states along the chain. We now see how this result connects with the previous analysis on the 
density relaxation times. In the small U/J regime the particle is able to escape from the well and to explore larger 
portions of the chain, thus resulting into a fast density relaxation. As opposite, for large [//J, the particle bounces 
back and forth inside the well, finding hard time to escape from it. This lack of diffusion results into a very long-time 
scale for the density to relax to its homogeneous value. 

The above results show explicitly that some kind of localization in the many-body configurational space does occur, 
at least in the finite system. 0, [Ij] While such an intriguing behavior might well be a subtle effect due to the finite 
size spectrum, it could also signal the onset of a genuine localization that survives in the thermodynamic limit. 

We conclude the discussion on the inliomogeneous initial states by studying the case at density n < 1, when at 
equihbrium there is no longer a Mott transition. For example, we consider n = 2/3 and an initial density profile 
(. . . 2, 0, 0, 2, 0, 0, . . . ). Interestingly, we find quite a different behavior for the density relaxation times th (sec Supple- 
mentary Material), with a much smoother crossover from small to large values of U/J and no evidence of any increase 
in the relaxation times with the system size. This fact suggests a deeper connection between the observed dynamical 
behavior and the zero temperature Mott transition that occurs at equilibrium and at integer filling, as suggested by 
calculations with infinite-coordination lattices. 



Homogeneous initial states and Quantum Quenches 

In light of the previous results, one may question that by choosing an inhomogeneous configuration of doublons 
we pick up a rather specific initial state in the Hilbert space. We are now going to show that the above findings 
strongly affect the dynamics starting from a perfectly homogeneous state. In this respect, a particularly interesting 
class of initial states are the ground states of H for given values of the interaction Ui, which are let evolve under 
the Hamiltonian dynamics after a sudden increase of the interaction to a final value Uf > Ui, the so-called quantum 
quench. KoUath and coworkers [2^ reported evidence for the existence of two separated regimes in which either 
thermal or non-thermal behavior is observed for local observables. The origin of the non-thermal behavior in the large 
Uf region and the possibility of an ergodicity breaking in the thermodynamic limit is still highly debated 2^ 3^. 
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In the following, we focus on an average density n = 1 and we show that signatures of long lived metastable 
states of doblons can be identified in the dynamics after a quantum quench. At variance with the previous numerical 
calculations, now both the initial state and the quantum Hamiltonian do preserve the spatial homogeneity and, 
therefore, the quest for possible signatures of ergodicity breaking requires a different approach. Since we have identified 
density relaxation as the slowest process in the problem, we monitor the dynamics of the system by measuring the 
auto-correlation of the density averaged over all sites, namely through 

C(^) = ^ E Mt)n^m {n,m . (3) 

i 

For any finite size system, limt_>ooC(t) = since the densities decorrelate at very long times. Indeed, as shown in 
Fig. m for small Uf <^ Ui this quantity has a very fast transient to zero. On the contrary, for Uf '3> Ui the density 
auto-correlation C{t) gets stuck into a long-lived finite value plateau ^ 0, before approaching zero only on a much 
longer time scale. If we extract a relaxation time from C{t), we find a similar behavior as in Fig. [2l i.e., a dramatic 
increase of the relaxation times above a threshold value of the final interaction strength. 

In agreement with the previous analysis, the appearance of such a long-lived metastable state characterized by a 
finite plateau of the density auto-correlation function might indicate an excess of double occupancies that have no 
channel to relax. In other words, the dynamical constraints brought by the interaction severely slow down density 
excitations, whose characteristic time scales increase abruptly after a critical threshold. The main phcnomenological 
traits of this dynamical arrest characterized by long-lived inhomogeneous states closely remind the physics of glassy 
materials. 



Variational description, lack of thermalization and higher dimensions 



The previous discussion revealed the existence of a threshold energy above which a steep increase of the relaxation 
time of density fluctuations takes place. In order to assess the relevance of this phenomenon for the dynamics of 
larger systems or even for higher spatial dimensions, it is desirable to devise a comprehensive alternative framework 
able to catch its very characteristics. Here, we introduce an approach based on real-time variational Monte Carlo 
(see Supplementary Material for details), which has two important advantages: it allows us to follow the evolution for 



times comparable to those accessible experimentally, which are much longer than t-DMRG; [31|, |321 it can be easily 
extended to higher dimensions. We mention that mean-field-likc variational approaches to the real-time dynamics of 
correlated systems have been developed in recent years. [13, [l^ 

However, although they seem to capture well the main features of the dynamical evolution, these methods are 
unable to describe important aspects such as damping and relaxation. On the contrary, our approach is sufficiently 
rich to account for damping and relaxation of local observables. It is based on a very simple and transparent out-of- 
equilibrium extension of the Jastrow-like variational wave function that was shown to describe quite accurately the 



Non-Ergodic 




Figure 4: Inverse relaxation times of density excitations in the homogeneous system, see Eq. ([3|. From left to right, different 
curves correspond to different initial states at Ui/J = 0, 1, and 2. Insets: real part of the density correlations C{t) in the 
ergodic and in the non-ergodic region. Data are obtained with exact-diagonalization on a lattice with A'^ = 12. 
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Figure 5: One-dimensional results for the time-dependent expectation values of the on-site potential energy £p{t) = 
-J- {ni(ni — 1)) in the ergodic {upper panel) and in the non-ergodic regions {lower panel). The initial state is the ground 
state of the Bose-Hubbard Hamiltonian with Ui — 2J and the considered system size is A'' = 200. Grand-canonical thermal 
averages are shown for comparison as dashed horizontal lines. 



equilibrium phase diagram of the Bose-Hubbard model: [33 



(4) 



where |5'o) is the initial state and Vij{ni, rij] t) is a Jastrow factor that depends on the occupancies n,; and rij of two 
sites i and j and varies with time so to maintain the time evolution as close as possible to the true evolution via the 
Schroedinger equation (see Supplementary Material) . The comparison between our approach and the t-DMRG SJ] is 
reported in the Supplementary Materials, demonstrating the high accuracy of the time-dependent variational Monte 
Carlo. 

Results for the time evolution of a local observable, such as the potential energy, after a sudden quench from Ui = 2 J 
to a final U f are shown in Fig. [S] The values of the thermal averages have been computed in the grand-canonical 
ensemble by means of finite-temperature quantum Monte Carlo calculations, (35| with the effective temperature fixed 
by the average energy of the initial state, which we take as the best variational approximation for the ground state at 
U, = 2 J. 

As shown in Fig.[5l in the region of small Uf we observe a damping of the average potential-energy, which approaches 
a quasi-steady stationary value in contrast to the simple Gutzwiller wave function, [2^ . [25j . In this regime, damping 
is mainly due to a density-density Jastrow factor of the form Vij{ni, fij]) = Vij{t) rii nj, which already at equilibrium 



was shown to provide a satisfactory description of the physical behavior. [331 This fact enlightens the relevance of the 
Bogoliubov modes whose dephasing during the time evolution allows to approach the stationary state. Remarkably, 
the steady state averages coincide with the thermal ones; a signal that the dynamics is ergodic. 

In the region of large interactions Uf, a simple density-density Jastrow factor does not account for all relaxation 
pathways, which will now mainly result from specific correlations among doublons, holons and between holons and 
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Figure 6: Two dimensional results for the time-dependent expectation values of the on-site potential energy £p{t) = 
{ni(ni — 1)) in the ergodic (upper panel) and in the non-ergodic regions (lower panel). The initial state is the ground 
state of the Bose-Hubbard Hamiltonian with Ui = 4J and the considered system size is A'^ = 20 x 20. Grand-canonical thermal 
averages are shown for comparison as dashed horizontal lines. 

doublons. The effective Hamiltonian ([2]) indeed explicitly shows that doublons attract each other as well as holons 
do, while doublons repel holons. These correlations, as well as other among higher on-site densities, can be easily 
implemented via the the Jastrow factor in Eq. (j4]) and indeed substantially improve the dynamics. Interestingly, the 
effective interaction between doblons that results from the dynamical variational calculation turns to be attractive, 
therefore leading to a consistent determination of the anticipated dynamical effects that drive the dynamics in this 
regime. As we see from Fig. [3 in the region of very large U f the potential-energy expectation values do show a 
damping to a non-thermal quasi-steady state on a time of the order of T£) ~ ^/Uf. This fast time scale must be put 
in comparison with the much longer one, of Fig. 21 It is natural to identify th with the time scale that controls 
the eventual escape from the quasi-steady state, hence the approach to thermal equilibrium. Whether this time scale 
does truly diverge in the thermodynamic limit, or rather saturates to a very large value which is still out of reach 
for state of the art numerics, it is certainly an important issue that cannot be definitively solved. However, we can 
safely state that a large finite system of actual experimental relevance will get stuck for a quite long time into highly 
inhomogeneous metastable states, which we revealed to be on the verge of a spatial symmetry breaking. 

The above results have been obtained for out-of-equilibrium one dimensional systems and therefore leave the ques- 
tions concerning the dependence on the dimensionality of the problem still open. However, the anomalously long-time 
relaxation of the density auto-correlation points towards a kind of glassy behavior that should be observable even in 
higher dimensions, provided that the interaction induces sufficiently strong dynamical constraints. In this regard, we 
have studied the two-dimensional case by means of our time-dependent variational scheme, verifying that a similar 
behavior occurs even in two dimensions. In Fig. [51 we show the results of the potential-energy expectation values as 
a function of the time. As before, if the repulsion is weak, the time average coincides with the thermal ones, while for 
strong repulsion there is a clear difference between thermal and time averages, giving rise to a nonergodic dynamics. 
We therefore argue that the phenomenology we have hereby identified is almost independent on the dimensionality 
and it is rather due to the existence in strongly correlated systems of high-energy incoherent excitations that do not 
have channels to relax efficiently. This scenario is also consistent with the observed strong doping dependence of the 
relaxation time (see Supplementary Material). 
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Discussion 

In conclusion, we have found that the dynamical constraints brought by a strong interaction can trap the evolution 
of repulsive bosons hopping on a lattice into metastable states that lack translational symmetry, provided that the 
energy stored into the initial state is above a threshold. We pointed out that a self-induced effective attraction among 
doublons is one of the major processes that can effectively freeze the dynamics on long time scales. Such a mechanism 
is recognized to play a role in the density relaxation processes of purely homogeneous systems through a dynamical 
arrest visible in time-dependent density correlations. The main features of this intriguing behavior, namely the slowing 
down of density excitations and the long-lived inhomogcncous pattern, resembles closely a kind of glass transition. 

Moreover, we have shown that the time evolution of the many-body problem can be mapped onto that of a particle 
moving from the edge of a semi-infinite tight-binding chain with nearest-neighbor hopping, where each site represents 
a many-body wave function. This model looks like an Anderson model, since both the on-site energy and the hopping 
vary from site to site, with a potential well at one edge due to the high-energy content of the initial state. Interestingly, 
we find a delocalization-localization transition in this problem, with the particle being unable to diffuse on the whole 
chain above a certain value of the well depth. We consider this analogy quite suggestive and potentially constitute an 
even stronger indication of ergodicity breaking in the many-body space which is worth to be further investigated. 
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—Supplementary Material- 
Exact Diagonalization 

The exact time evolution of an initial quantum state |^'(0)) is given by: 

|*(t)) = e-^« *|*(0)) = 'l'^'")' (5) 

n 

where |(f>„) are the eigenstates of H, with corresponding eigenvalues £"„, and c„ = ($„|5'(0)). Therefore, in order 
to follow the quantum dynamics to large times, the full knowledge of the spectrum is in principle needed. However, 
whenever the Hilbert space is so large that it does not allow for a full diagonalization, an alternative approach is in 
order. At each time, we can make an evolution for a time step At by considering a truncated Taylor expansion of the 
unitary evolution: 

|*(t + Ai)) =e^'«^*|*(t)) ~ ^ lz!^7^fc|^(t)). (6) 

Here, the series converges very quickly with k and the cutoff kmax must be chosen to obtain the desirable convergence 
for a given At. Therefore, the evolved wave function + At)) can be easily obtained by summing terms that can 
be in turn recovered by repeatedly applying the Hamiltonian H to |\E'(t)). The full time evolution at long times can 
be achieved by subsequent small-time evolutions. 

Moreover, time-dependent correlation functions of the form CA{t) = (^|A^^(t)A(0)|^'), where A is an arbitrary 
operator, can be as well obtained by first applying A to then performing the time evolution \'^A{t)) = e^'^ '^l^*), 
and finally computing C^(t) = {<ij{t)\A''\^Ait)), where |*(t)) = e"*^ *|*). 

The numerical accuracy of the unitary evolution can be verified by checking the conserved quantities of the time- 
evolution (for example the total energy) which in our calculations remain to all purposes constant up to the longest 
considered evolution times. 



Lanczos Method 



In this work we have used the Lanczos method to construct a basis of L orthonormal states in which the full 
Hamiltonian takes a reduced tridiagonal form: 



/ 



-Hi 



£0 ^0 







t2^1 £2 



V 









(7) 



The basis set is constructed in a recursive way starting from an initial vector |0) upon repeatedly applying the 
Hamiltonian 



tn^n+i \n + 1) = U\n) - e„ \n) - t^-i^n \n - 1) , 

with the coefficients of the tridiagonal matrix given by 

tn^n+i = (n+l|H|n) 
e„ = {n\H\n). 



(8) 
(9) 



Once the tridiagonal matrix is constructed we are in condition to consider a time-dependent problem of a particle 
initially localized on site |0) whose dynamics is governed hy Hl- 

A final important remark concerns the practical numerical evaluation of the Lanczos states and matrix elements on 
calculators with a finite-precision arithmetic. Numerical truncation errors can indeed accumulate and lead the most 
recently generated Lanczos vector not to be orthogonal to the previously generated ones. To overcome this problem, 
we have carefully analyzed the dependence of the results on the machine precision, considering up to 10'* bits of 
floating point arithmetic, for which no significant truncation error manifests up to the largest value of L considered 
in the paper. 
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Time-Evolving Block decimation 



The general idea of the Time-Evolving Block decimation algorithm |15l . |17| is to consider a particular representation 
of a generic many-body state defined in a Hilbert space of dimension , where V is the dimension of a certain local- 
basis. Such a generic state is fully described in terms of the local quantum numbers ik as 



V 

1=1 

and in the Vidal's representation the coefficients Ci^i2..iN ^'^^ taken to be 

X 

C- . . = V AWrW^i A^rl^l^^ \[3]p[3]i3 \[4] . _y\N] p[Ar]iiv \[N+l] qqn 

where the local tensors F and the bond matrices A constitute a set of Vy^N + x + 1) parameters that specify 
the state in this representation. The time evolution of a given initial state is obtained upon considering a small 
time-step At, and a repeated application of the real-time propagator e~'^^* within a Suzuki- Trotter approximation 
scheme, inducing systematic errors of magnitude O (Ai"^). When the infinitesimal unitary evolution is applied the 
decomposition (jlOp is consistently updated. The other major source of systematic errors is however intrinsically due 
to the finite amount of entanglement retained by Vidal's representation. We have checked the convergence of the 
studied properties with respect to both the time step and the entanglement cutoff parameter, considering a maximum 
of X = 400, for which the density is converged up to the plotted times of Fig. 2 (right panel) in the main paper. 



Real-time variational Monte Carlo 



In this part of the Supplementary Materials, wc give the details of the real-time variational Monte Carlo approach. 
First of all, we would like to stress that this method docs not suffer from intrinsic instabilities arising from the sign (or 
phase) problem. Indeed, our method is essentially a straightforward generalization of the variational technique, where 
the probability density is strictly non-negative. Although approximated, this approach has firm physical basis and 
can accurately reproduce numerically exact results in cases where t-DMRG method can be applied (e.g., mostly in 
one-dimensional lattices). Nonetheless, variational Monte Carlo can be generalized and applied in larger dimensions 
without any further numerical instability. 

In general, it is reasonable to expect that the relevant physical properties driving the dynamics can be reduced to 
a number of parameters significantly smaller than the size of the full Hilbert space of the problem, as much as it is 
done in renormalization-based numerical and analytical tools. Let us introduce a set of "excitation operators" Ok, 
which are expected to describe the relevant excitations over the ground state that are involved in the dynamics of the 
system. The actual set of these operators has to be carefully chosen, depending upon the properties of the quantum 
phase. In the case of one-dimensional interacting bosons. Ok can be selected to be density-density correlations, e.g.. 
Ok = 'Y^^niiiii+k^ or doublon-doublon, holon-holon, and doublon-holon correlations. As discussed in the text, the 
presence of the former ones is sufficient to obtain a very accurate time dynamics for small quenches, whereas the 
inclusion of the latter ones is necessary for large quenches. In Fig. [3 we report the comparison of the t-DMRG 
dynamics with our variational Monte Carlo approach. 

Let us discuss the details of our method. First of all, we define a basis |x) in which the occupation numbers of 
the bosons are defined in the lattice. Since the operators Ok arc diagonal in this basis, we will define Ofe(x) as the 
(time-independent) value that the operator takes over the configuration |x). The variational wave function is given 
by: 



*(x,t) = e^(*^exp 



Y,Ok{^)ak(t) 



e(x), (11) 



where £,{x) is a time-independent uncorrelated state, such as the bosonic condensate, whereas 9{t) and ak{t) are 
time-dependent complex parameters. 

The variational form introduced here allows us to map the quantum dynamics associated to a given Hamiltonian 
% into the dynamics of the variational parameters 9{t) and ak{t). The operators Ok can be thought as collective 
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Figure 7: Time-dependent expectation values of the on-site potential energy £p{t) = ^ {ni{ni — 1)). The initial state is the 
ground state of the Bose-Hubbard Hamiltonian with Ui = 2J. Real-time Variational Monte Carlo results are obtained for an 
N = 100 chain with periodic boundaries. Corresponding data by t-DMRG [sj are the potential energy at the center of an 
N = 64: chain with open boundaries. 



variables for the quantum dynamics; on physical grounds, we expect that the quantum dynamics will be characterized 
by a set of collective variables whose number Np is much smaller than the original size of the Hilbert space. This 
reduction assumption is the key requirement for having an efficient time-dependent variational approach. 

The equations of motion for the variational parameters can be found by minimizing the Euclidean distance in 
the Hilbert space between the exact and the variational time dynamics. Given at a certain time the quantum state 
^{x,t), its exact time evolution leads to <j){x,t) = Vl'(x, i) x [— if(x, <)], where the complex-valued "local-energy" 
• other hand, an infinitesimal change in time in the variational parameters leads to 



is f(x, 

i'ix,t) = *(x,t)x e{t) 



J2k^k{^)ctk{t) ■ The minimization of their mutual distance = "/"(^i ~ ^(^i 



gives a set of coupled differential equations for dk{t) and 0{t): 

k' 



(12) 
(13) 



Ok 



- {Ok) and 

|^(x,f)|^F(x) 



denote expectation values over the square modulus of the wave function at time t, 



where 60k 

namely (F) — ^ , ,,,2 

The coupled equations of motion determine a consistent Hamiltonian dynamics in which at all times both the norm 
N{t) = |V'(x, t) ^ and the expectation value of the energy Eit) ~ {£{t)) are strictly conserved. Indeed, one can 

show that N(t) oc + Efc^fc (i)(Cfc)j and E{t) oc J2k {o^ki^OkS") + ai{SOk£^)}, which are both vanishing for 
the optimal solutions for which 2?(t) is at a minimum. 

Not surprisingly, the same equation of motions that we have derived by means of an optimal distance principle can 
be also derived from the variational principle, namely it can been shown that the action S = J dt{'i' ~ TL\'^) is 
indeed stationary for the time-trajectories of the variational parameters induced by Eqs. and p^ . 

The time-dependent expectation values over the square modulus of the many-body wave function cannot be calcu- 
lated analytically for generic operators Ok, and a variational quantum Monte Carlo is needed to evaluate them. In 
particular, |t/'(x, t)|^ can be straightforwardly interpreted as a probability distribution over the Hilbert space spanned 
by the configurations x and a Markov process can be devised, such that the stationary equilibrium distribution 
coincides with the desired probability measure. 



Inhomogeneous states dynamics at n — 2/3 



In the main paper, we have concentrated our attention on the unitary-filling density, where at equilibrium a Mott 
transition takes place upon increasing the interaction strength. To better realize the importance of the density on 
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Figure 8: Exact time evolution of the densities on initially doubly occupied sites at filling n — 2/3 for an = 15 chain. Different 
curves from left to right correspond to increasing values of the final interaction ai U/ J — 0,1, . . . 10. Inset: corresponding inverse 
relaxation times for increasing lattice size N = 9, 12, and 15. 

our considerations on the largc-C/ relaxation times, we now show the exact time evolution of inhomogeneous states at 
n < 1. For example, we consider n ~ 2/3 and an initial density profile (. . . 2, 0, 0, 2, 0, 0, . . . ). As shown in Fig. [SI we 
find quite a different behavior for the density relaxation times with respect to the unitary-filling case, with a much 
smoother crossover from small to large values oiUf. Moreover, there is no evidence of any increase in the relaxation 
times with the system size. This fact suggests that the dynamical constraints brought by the effective interaction 
among doblons are much stronger at n = 1, where a sharp crossover in the relaxation times is observed. We finally 
remark that, even at non unitary fillings, we expect larger clusters of doblons to have larger relaxation times than the 
simple initial state we have considered here. The possibility of a non-thermal behavior at non-unitary filling cannot 
be excluded, even though the strongest manifestation of the dynamical arrest is expected to happen at n = 1. 

Lanczos-basis analysis of an integrable model 

To better elucidate the many-body localization phenomenon suggested by the Lanczos-basis analysis of the Bose- 
Hubbard model, we hereby show that localization is not strictly implied by integrability. To this purpose, we consider 
a model of one-dimensional hard-core bosons prepared in an initial inhomogeneous state with an alternating density 
profile (...1,0,1,0,...) and let it evolve in the Lanczos basis with their non- interacting (integrable) kinetic Hamil- 
tonian. In Fig. ^ we show both the Lanczos hopping elements and the expectation value of the Lanczos particle 
traveling in the many-body space. The effect of integrability is to reduce the number of allowed Lanczos states with 
respect to the total Hilbert space. Indeed, we find that the hopping is a strongly decreasing function of the iterations, 
approaching zero after a certain number of states v{N) that is a finite fraction of the full Hilbert space. Noticeably, 
the average position traveled by the particle increases accordingly, leading, therefore, to a full delocalization in the 
thermodynamic limit. This analysis shows that a localization in the Lanczos basis is not necessarily due to integrabil- 
ity - whose effect only amounts to reduce the number of active many-body states - whereas it is due to the effective 
dynamical constraints brought by the interaction. 
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Figure 9: Left panel - Nearest neighbor hopping of the effective chain that represents the non-interacting hard-core bosons 
Hamiltonian in the Lanczos basis starting from (. . . 1, 0, 1, 0, . . . ) Right panel ~ Time-dependent expectation value of the wave- 
packet position traveling in the Hilbert space generated by a chain of viN) Lanczos states. The horizontal lines correspond to 
v{N)/2 and show that delocalization takes place. 



